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This paper presents a numerical analysis of the transition from selective withdrawal to 
viscous cntrainmcnt. In our model problem, an interface between two immiscible layers 
of equal viscosity is deformed by an axisymmetric withdrawal flow, which is driven by 
a point sink located some distance above the interface in the upper layer. We find that 
steady-state hump solutions, corresponding to selective withdrawal of liquid from the 
upper layer, cease to exist above a threshold withdrawal flux, and that this transition 
corresponds to a saddle-node bifurcation for the hump solutions. Numerical results on the 
shape evolution of the steady-state interface are compared against previous experimental 
measurements. We find good agreement where the data overlap. However, the numerical 
results' larger dynamic range allows us to show that the large increase in the curvature 
of the hump tip near transition is not consistent with an approach towards a power-law 

cusp shape, an interpretation previously suggested from inspection of the experimental 
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measurements alone. Instead the large increase in the curvature at the hump tip reflects 
a logarithmic coupling between the overall height of the hump and the curvature at the 
tip of the hump. 



1. Introduction 

Topological transitions of a fluid interface are a key step in many important and 
complex physical processes. Simple examples include the formation and coalescence of 
droplets where local stresses typically determine the fluid flows near the transition. Here, 
we examine a different scenario where the topological transition in the interface shape is 
driven by a large-scale steady-state flow. Figure [1] depicts the experiment which inspired 
our numerical study. A deep layer of viscous silic one oil overlays a seco nd, immiscible 



layer, comprised of a mixture of water and glycerin (jCohen fc Nag el (2002)). Steady-state 
large-scale flows are created in the two fluid layers by withdrawing liquid through a tube 
placed in the upper layer and replenishing the same amount of liquid back into the upper 
layer at locations far away from the point of withdrawal. The rate of fluid withdrawal Q 
dramatically affects the shape of the steady-state interface. When Q is below a threshold 
value Q c , the steady-state interface forms a hump. The lower-layer flow has the form of 
a toroidal recirculation, with a stagnation point at the tip of the hump. Since only liquid 
in the upper layer is withdrawn into the tube, this is known as the selective withdrawal 
regime. When the withdrawal flux Q is larger than Q Cl the flows and the interface remain 
steady, but the interface attains a different topology. The steady-state interface is now 
"open" and has the shape of a spout which extends smoothly into the tube. Since liquid 



is now withdrawn from both layers, this is known as the viscous entrainmcnt regime. In 
this regime, the amount of liquid in the lower layer slowly decreases while the amount 




Figure 1. (a) Schematic of experiment. Liquid is withdrawn at a prescribed volume flux Q 
through a tube inserted into the top layer and replenished at the same rate. The top layer is a 
viscous silicone oil. The bottom layer is a mixture of water and glycerin. The two layers have 
comparable viscosities. The flow rates are small so that viscous effects are significant in both 
layers. The tube diameter is 1.6 mm. (b) Selective withdrawal regime: below Q c a steady-state 
hump forms on the interface, (c) Viscous entrainment regime: above Q c a steady-state spout 
forms. Photos courtesy of Cohen & Nagel. 



of liquid in the upper layer slowly increases. The large-scale toroidal recirculation in the 
lower-layer is weakly perturbed by the onset of entrainment. The stagnation point is 
deflected slightly downwards, from the hump tip to the interior. 

Analogous transitions from selective withdrawal to entrainment occur in many indus- 
trial and physical processes. Some examples on large lengthscales incl ude the intrusion of 



water i nto an oil reservoir dur i ng th e last stages of petroleum recovery ( Renardv fc Joseph 



( 19851) : 



Renardv fc Renardvl (|l985l )). water-coning in sandy beds (jForbes et al 



the for mation of long-lived plumes by thermal convecti on in the earth's mantl e 



(|2002h ) and pollution recovery on lakes and oceans (jlmberger fc Hamblin 



small lengths c ales, the bursting of a liqu i d drop i n a st eady straining flow 



Grace! (jl9821) 



Rallison fc Acrivosl ([1978T I 



Stone 



tendril by a drop sliding down an incline 



(2004)). 



Jell nek fc Manga 



Tavloi 



1982)). On 



(| 19341 ) 



19941) h the depo sition of a thin trailing 



Limat fc Stond |2004J)) and the formation of 



thin cylindrical tendrils and droplets in microfluidic devices using axisymmctric now- 
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focusi ng (jAnna et al 



(2003); 



Link et al. 



(2004): 



Lorenceau et al. 



(2005); 



Utada et al 



(|2005l )) all involve topological transition of steay-state interfaces. 



In the selective withdrawal experiments, the cur vature at the hump tip i ncrease s signi f- 



Cohen fe Nagel 



420 02J) 



Cohenl (120041) . 



icantly as the transition is approached from below 
Intriguingly, measurements of the hump height and curvature evolution near the transi- 
tion suggest the steady-state hump shape evolves towards a power-law cusp. In order to 
explain this apparent divergence, Cohen and Nagel made an analogy to recent works on 
topological transitions that correspond to the formation of a finite-time singularity in the 



govern ing equations. S ome example s are t h e break-up of a l i quid d r op into two 



drops (jEggersI (jl997l ); 



Cohen et al 



U2m 



the coalescence of two d rops into one 



Thoroddsen et al 



Zhang fc Lister 



Qguz fe Prosperetti 



1999) 



19901) 



d aught er 



Doshi et al 



Eggers et al. 



(2003)), 



19991) : 



(|2005f )). and the eruptio n of an electrohydrodynami c spout under the 



sudden application of a large electric- field (jOddershede fc Nagell (|2000| )). In all these sit 



uations, the topology change necessarily requires a divergence in the mean curvature at 
a point on the interface. As a result, the Laplace pressure due to surface tension, which is 
proportional to the mean curvature, diverges as well, thereby creating a singularity in the 
governing equations. This suggested analogy between the topological transition in vis- 
cous withdrawal and finite-time singularity formation is surprising. While the break-up 
of a liquid drop requires that the interface deforms continuously from a single, connected 
shape to several disconnected shapes, and therefore requires the formation of a singular 
shape at break-up, the topological transition in this experiment deals with changes in 
the steady-state interface shape. Therefore the transition from selective withdrawal to 
viscous entrainment does not require that the hump shape evolve into the spout shape 
continuously as a function of the withdrawal flux Q. Generically, one would expect a 
discontinuous change in the steady-state shape across the selective withdrawal to viscous 



cntrainment threshold. The apparently, nearly continous shape-change observed in the 
experiment is, therefore, rather unusual. 

This surprise motivates our work. Here we conduct a numerical analysis of a model 
withdrawal process in which the interface deformation is solely controlled by viscous 
stresses and surface tension. We focus on the selective withdrawal regime and analyze 
how the interface fails at the transition. In particular we arc interested in elucidating 
the various factors that determine the final hump shape. Experiments show that this 
final shape is well correlated with the minimum spout thickness attained in the viscous 
entrainment regime (figure [TJ;) . Consequently, a detailed understanding of this topolog- 
ical transition will help advance a variety of practical applications. One such applica- 
tion entails takeing advantage of viscou s entrainment to encap sulate biological cells for 



transplant therapy (jCohen et al 



(120011) : 



Wvman et al 



(|200J)). In this technique, the 



minimum spout size directly determines the minimum thickness of a protective polymer 

coating that can be applied to the cells. 

The rest of this paper is organized as follows: In section 2 we discuss related works 

and background materials. In section 3 we describe the model withdrawal process and 

its numerical formulation for a system in which the upper and lower fluid viscosities are 

equal. In section 4 we analyse solutions to the model. Section 4.1 describes in detail the 

interface evolution near the transition and the sharpening of the hump tip. We find that 

scalings of the hump curvature and hump height with the rate of withdrawal identify 

the transition as a saddle-node bifurcation. In accordance with this picture we find that 

the hump solution remains smooth as the transition is approached. Surprisingly, we also 

find that near the transition, the hump height couples logarithmically to the curvature 

at the hump tip. In section 4.2, we show that, for the equal viscosity layers investigated, 

these key findings are independent of the particular withdrawal conditions in the numer- 
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ical simulations. In section 4.3, we compare the numerical and experimental results and 
find excellent agreement between them. Moreover, we show that the unusual logarithmic 
coupling between the hump height and the hump curvature can easily lead to a misin- 
terpretation that the system approaches a steady state singularity at the transition. A 
more detailed discussion of the distinction between the previous interpretation and our 
results is presented in section 5. 

Our findings indicate the transition from selective withdrawal to viscous entrainment 
involves a discontinuous change in the steady-state shape that has the structure of a 
saddle-node bifurcation. This discontinuity is effectively obscured in practice by the loga- 
rithmic coupling between the hump height and the hump curvature, which causes changes 
in the hump curvature to be far larger, and therefore more evident, than changes in the 
hump height as the transition is approached. 



2. Background 

Many different realizations of flow-driven topological transitions have been investigated 
in previous works. Here, we review those studies that directly assess the ingredients neces- 
sary for a topological transition to take place and the appropriate ways of characterizing 
the evolution of the interface near such a transition. 

2.1. Viscous entrainment in the absence of surface tension 

A number of experimental and theoretical works have focused on the transition from 
selective withdrawal to viscous entrainment in stratified, two-layer systems where surface 
tension effects are assumed to be either negligible or weak. These studies were motivated 
by geophysical flows, in particular the emptying of a magma chamber during a volcanic 
explosion. Blake & Ivey investigated the threshold withdrawal fl u x Q c necessary for the 



onset of entrainment in miscible, stratified layers (jlvev fc Blakd (|1985[ )). In a follow up 



numerical study (jListerl (|1989l )). Lister analyzed the onset of entrainment in two liquid 
layers of equal viscosity and found that a finite value for Q c exists only if surface tension 
effects, however small, are included in the analysis. These studies did not analyze the 
hump evolution near the transition. 



2.2. 2-D analogues 



Air entrainment, such as occcurs in high-speed coating (jSimpkins fc Kuck ( 200C )) and in 



the im pact of a jet of viscous liquid into a layer of the same viscous liquid 



Lorenceau et al 



( 20041 )) is o ften accompanied by the interface approaching a two-dimensional power-law 



cusp shape ijJoseph et al 



( 19911 )). Jeong & Moffatt showed analytically that, in the ide- 



alized limit where air flow eff ects are negligibl e , ther e is no transition from selective 



withdrawal to air entrainment (jJeong fc Mo ffatt 



19921 )). Instead the interface shape ap- 



proaches a steady-state singularity in the shape of a power-law cusp as the fl ow rate 



is incr e ased towards infinity. M ore recent thereotical and experimental analyses ([Eggers 



(200l|); 



Lorenceau et 



(200,3)) showed that, when air flow effects are included, they per- 
turb the idealized dynamics so that a transition to air entrainment occurs at finite flow 
rate. The approach towards a steady-state singularity is then cut-off at a small length- 
scale. Since the cut-off for air entrainment originates from viscous stresses associated 
with air flow, the cut-off lengthscale has a strong dependence on the viscosity contrast 
and increases as (/Uq/a*) 4 ^ 3 , where /iq is the air viscosity. 



2.3. Selective withdrawal in 3D with viscous flow and surface tension 



In the ir studies of 3-D axisymetric selective withdrawal, Cohen & Nagel (jCohen fc Nagel 



((2002J)) used two liquids of comparable viscosities and measured h, the hump height, 



and k, the mean curvature at the hump t ip, as a function of the withdrawal flux Q. 



In a follow up study Cohen (jCohenl (|2004f )) analyzed the transition for different pairs 
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of liquids. In particular, the viscosity ratio, defined as the viscosity of the liquid to be 
entrained divided by the viscosity of the entraining liquid, was varied from 0(1) to 10 -3 . 
Analysis of the measurements show that, as the transition from selective withdrawal to 
viscous entrainment is approached from below, the steady state hump curvature increases 
dramatically. Cohen and Nagel showed that this rate of increase can be fit to a power 
law divergence. They also noted however, that the hump never reaches a cusp shape. 
Instead the final stage of the hump evolution is cut off by the transition. This behavior 
lead to the interpretation that a singular solution for the steady-state shape organizes 
the nearly continuous transition between the steady-state hump and spout shapes. The 
transition from a hump to a spout was observed to occur when the radius of curvature 
at the hump tip was 0(50) fxm or more. Moreover, the cut-off lengthscale shows little 
dependence on the viscosity contrast between the two layers. Consequently axisymmctric 
viscous withdrawal in 3D diffcrcs strikingly from air entrainment in 2D, which shows a 
strong dependence of the cut-off lengthscale on the viscosity contrast. 

2.4. Viscous drainage 

An analogous, but not necessarily equivalent, topology change occurs in viscous drainage 
experiments, in which a viscous liquid exposed to air is placed in a container with an 
opening at the bottom surface and is allowed to drain out o f the container due to its own 



weight |Chaiebl |2004j); 



Courrech du Pont fc Eggersl (|2006l )) . As in the setup by Cohen 



& Nagel, the drainage creates a large-scale toroidal recirculation whose viscous stresses 
deform the interface. Thus the drainage experiment resemble an upside-down version 
of the two-layer withdrawl experiment. A key difference is that the layer depth of the 
withdrawn fluid changes with time in the drainage experiment, but remains constant in 
the setup used by Cohen & Nagel. As the liquid drains out and the layer depth decreases 
below a critical value, a sharp cusp develops on the interface, a feature identified by 



Correch du Pont & Eggers as a steady-state singularity. A further decrease in the layer 
depth causes the cusp to approach closer to the bottom surface and eventually intrude 
into the aperture through which the fluid is draining. Correch du Pont & Eggers report 
that there is no transition from selective withdrawal to viscous entrainment throughout 
the entire drainage process. Moreover, the maximum curvature obtained in the drainage 
experiments, which corresponds to using two fluid layers with a viscosity contrast of 
10 -6 , are far larger than the maximum curvature measured by Cohen & Nagel for two 
fluid layers with a viscosity contrast of 10~ 3 . The origin of these different outcomes is, 
at present, not understood. 

2.5. Viscous withdrawal in 3D with surface tension 

Motivated by the thin stable spout created in the viscous entrainment regime, Zhang 
analyzed the transition in the reverse direction, from viscous entrainment to selective 
withdrawal. The analysis focuses on the regime where the viscosity of the fluid entrained 
is far smaller than the viscosity of the ent r aining liquid, and constrains the spout shape to 



be nearly cylindrical throughout (|Zhand (|2004l )). These simplifications allow the steady- 



state spout shape to be described via a long-wavelength model. Surprisingly, results 
for the model show that, changing the boundary conditions on the interface, without 
changing the material parameters, can profoundly change the nature of the shape tran- 
sition. For some boundary conditions, the steady-state interface deforms continuously 
with withdrawal flux from a spout to a hump across the transition. For other boundary 
conditions, the steady-state interface changes discontinuously at the transition. Thus, 
there appears to be a subtle interplay between the entraining flow and surface tension 
effects. As a result, the transition is strongly influenced by constraints on the large-scale 
shape of the interface. Case fc Nagel have recently characterized the steady-state spout 



shape experimentally (jCase fc Nagell (|2006l )^ 



2.6. Emulsification 

Another extensively studied shape transition that has strong points of similarity with 
the selective withdrawal to viscous entrainment transition occurs when a singl e liqui d 



drop is emulsified into se v eral drops by an axisymmet ric straining flow (| Taylor] (|1934T ); 



Rallison fc Acrivod (| 19781 ): 



Grace! rtl982l) 



Stond |l99J)). In both situations, an imposed 



flow far from the interface induces flows near the interface between two viscous liquids. 
The steady-state interface is deformed at low flow rates and "broken" at high flow rates. 
The crucial differences are that, for drop emulsification, there is no steady-state interface 
above the threshold flow rate. The drop is simply stretched out indefinitely by the flow. 
Also, even in the analog of the selective withdrawal regime, when a steady-state shape 
exists for the drop, the flow dynamics are different because the flow inside the drop 
must assume a form satisfying global volume conservation. In contrast, for the two-layer 
withdrawal experiment, the layer depth is always much larger than the hump height. 
It is then reasonable to expect that the global volume constraint does not introduce a 
significant modification to the flow dynamics. We will show this expectation is indeed 
borne out by results from the numerics. 

The viscosity contrast between the drop and surrounding fluid strongly affects the 
deformation-to-burst transition. When the drop is far less viscous than the surrounding 
liquid, the steady-state shape approaches a cusp shape as the burst transition is ap- 
proached from below. In particular, the radius of curvature at t he two, elonga t ed en ds 



of the liquid drop is cut-off on an exponentially small lengthscale (jAcrivos fe Lol (|1978l )). 
Long-wavelength analyses of the extended drop shape in this regime also show that the 
burst transition correspond s to a saddle- n ode bifurcation in the leading-order solution 



for the overall drop shape (|Tavlorl (|19641 ): 



Buckmasterl (|1973l )). When the drop viscos- 



ity is comparable with the surrounding liquid, the steady-state shape remains nearly 
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spherical as the burst transition is approached. Recently, numerics and phcnomcnologi- 
cal arguments have shown that the shape transition when the drop visco sity is equal t o 



the surrounding liqu i d visc osity corresponds to a saddle-node bifurcation (jNavot 



Blawzdziewicz et al 



( 20021) ) 



19991) : 



2.7. Inviscid selective withdrawal 



Finally, we note that a similar flow-driven topologi cal transition exis t s at high Reynolds 



numb er when water drains out of a filled bathtub (jSautreauxl (|1901l) : 



Lubin fc Springer 



(|l967T )). In practice, this is always accompanied by the formation of a vortex. More gen- 
erally, withdrawal from two stratified layers of inviscid fluids is relevant for transport 
in wat er storage reservoirs, where often a layer of fresh water overlies a layer of salty 



water (jlmberger fc Hamblinl ((1982)). Intriguingly, some of the results for the idealized 



problem of withdrawal from two l a yers o f perfectly inviscid fluids ( Tuck fc Vanden-Broeck 



(119841) : 



Vanden-Broeck fc Keller 



19871 ): 



Miloh fc Tvvandl (| 19931 )) suggest that, for the 



inviscid fluid system, the transition from selective withdrawal to full entrainment corre- 
sponds to the formation of a singular shape on the steady-state interface. This has moti- 
vated many recent theoretical and numerical studies, despite the fundamental difficulty 
that, without the small-scale smoothing provided by viscosity and/or surface tension in 
real life, the idealized problem is fundamentally ill-posed and can only be approached via 
careful limiting processes. For the same reason, direct comparison with experiments has 
also been difficult. For our study, the most relevant and suggestive results are that two 



dimen si onal withdrawal diff ers significantly from axisymmetric withdrawal (jForbes et al 



d200J) 



Stokes et al 



(|2005l )). Also, in both 2D and axisymmetric withdrawal, whether 
selective withdrawal or full entrainment is realized can depend crucially on the initial 
conditions which determine transient ev olution of the interface, i nstead of being solely 



determined by the boundary conditions (jHocking fc Forbea (|2001l )) 



In sum this review of viscous flow driven topological transitions shows studies of flow- 
driven topological transitions have produced a wealth of surprises. The wide range of 
behaviors observed indicate that apparently trivial differences in the exact experimen- 
tal configurations can have a dramatic effect on the transition structure. On the other 
hand, the interface evolution observed in these different situations shares a number of 
common features, such as the existence of a topological transition at a finite flow-rate, 
and the development of small-scale features on the interface near the transition. In order 
to make progress on the 3-D selective withdrawal problem, it is likely that experiments, 
simulations and theoretical calculations, all of which address the same flow geometry, 
will be needed. This is one reason why our study focuses on devising and analyzing 
a mod el withdrawal process that can be directly compared with the previous experi- 



ments (|Cohen fc Nagell (|2002l )) 



3. Modeling Selective Withdrawal 

This section describes our numerical model in the following order: we first identify 
the key parameters in the viscous withdrawal experiments by Cohen & Nagel and then 
describe an idealized withdrawal problem which retains these key features. A numeri- 
cal formulation of the idealized problem, together with the approximations, is given in 
section 3.1. Section 3.2 gives details of the numerical implementation. The key approxi- 
mations are that in our calculation the gradual flattening of the liquid interface due to 
hydrostatic pressure on large lengthscales is approximated by a hard-cutoff lcngthscale 
a. The gradual decay of the induced flow in the lower layer below the interface is ap- 
proximated by a cut-off condition requiring that the pressure in the lower layer becomes 
uniform laterally across the liquid layer. Comparisons of results obtained using these ap- 
proximations against those calculated using more realistic boundary conditions (section 
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4.2) and the experiment (section 4.3) show that the key features of the withdrawal are 
well reproduced by the simple numerical model. 

In the experiment, the height of the hump created is at most a few millimeters. The 
tube height S p is comparable in size, ranging from 0.255 cm to 0.830 cm. The capillary 
lengthscale £ 7 , which characterizes the relative importance of surface tension and hy- 



drostatic pressure, is about 0.55 cm. Deformations on lcngthscales below £ 7 = iJj/Apg 
arc stabilized by surface tension. Deformations on lengthscalcs above £j are stabilized 
by hydrostatic pressure. All these lcngthscales are much smaller than the dimensions of 
the liquid layer which is 12 cm deep and 30 cm in extent. Also, when tube height is 
sufficiently large, the interface deformation changes very little when the tube diameter 
is changed. This is the regime in which all the experimental data were taken. 

Since the tube height S p and the capillary lengthscale £ 1 are comparable, both length- 
scales are expected to influence the dynamics. As a result, there are several reasonable 
ways to define the Reynolds number, and it is unclear how the Reynolds number of the 
withdrawal flow should change with S p or £ 7 . To avoid this ambiguity, we use the ob- 
served hump height h as the characteristic lengthscale and define the Reynolds number 
Re as pQh/([pLTrSp), where p is the density of the upper liquid, p the viscosity of the 
upper layer, and Q/4wS p a typical flow rate. An estimate using parameters from the 
experiment yields Re m 0.1, indicating that incrtial effects arc small. The stress balance 
on the interface is therefore characterized by the capillary number, 



M Q 
7 4tt52 



Ca=^ (3.1) 



a ratio of the strength of the viscous stresses exerted by the flow relative to the Laplace 
pressure due to surface tension. For a typical experiment, the transition from selective 
withdrawal to viscous entrainment corresponds to a threshold capillary number Ca c in 
the range of 10 -3 . Thus surface tension effects are strong even at the transition. 



3.1. Minimal Model 

These results from the two-layer viscous withdrawal experiment motivate the following 
idealization. Since the container size is much larger than the characteristic lengthscale of 
the hump, the two immiscible liquid layers are taken as infinite in extent and in depth. 
We use a cylindrical coordinate system where z = corresponds to the height of the 
undisturbed, flat interface and r = is the centerline of the withdrawal flow. To drive 
the withdrawal, we prescribe a sink flow 

Uext(x) = i ^(xs - x) (3.2) 

where x s = Se z corresponds to a sink placed at height S above the undisturbed interface 
and Q is the strength of the point sink. This choice allows S to be the only lengthscale 
imposed on the flow, consistent with the experimental observation that the hump shape 
depends primarily on the height of the withdrawal tube S p and not on its diameter. 
Instead of explicitly including the effect of hydrostatic pressure, we require that the 
interface is deformed by the withdrawal flow only within the region < r < a. The 
interface Si is pinned at deflection at r — a. Beyond this point, the interface is taken 
to be entirely flat. Physically this hard cut-off, or pinning length a corresponds roughly 
to the capillary lengthscale i 1 . In the calculation the slope of of the interface at r = a 
is adjusted so that the pinning condition is satisfied at all flow rates. This introduces an 
error in the interface shape near r = a. This error is guaranteed to be small when the sink 
height S is far smaller than a. It turns out to be also small even when S is comparable 
with a. (In the experiment, tube height S p is cither comparable with or smaller than the 
capillary lengthscale £ 7 .) 

Finally, since the Reynolds number associated with flows in the experiment is small 
and the interface evolution is observed to remain essentially the same even as the layer 
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viscosities are made unequal, we idealize the flow dynamics in both liquids to be purely 
viscous. We also assume the two liquids have exactly the same viscosity. These assump- 
tions mean that the flow dynamics in the upper layer satisfies 



V ■ en = -Wpi +^V 2 Ui = 0, V-ui = 



(3.3) 



where Ui is the disturbance velocity field created in the upper layer due to the presence 
of the liquid interface, pi is the pressure and cri is the stress field associated with the 
disturbance velocity field Ux- Note u ext is a solution of Laplace's equation and therefore 
does not give rise to a uniform pressure distrubtion. The flow in the lower layer is created 
solely via the flow in the upper layer past the interface. Therefore, U2, the disturbance 
velocity field in the lower layer, satisfies 



V ■ er 2 = -Wp 2 + mV 2 u 2 = 0, V ■ u 2 = 



(3.4) 



where p2 is the pressure field in the lower layer, cr 2 the fluid stress field in the lower layer. 

Since the flow dynamics in both layers are completely viscous, and therefore described 
by the linear Stokes equations, we can represent the velocity field as an integral over a 
suita bly defin e d clos e d surface, as dis c ussed i n earlier works an d in textbooks on viscous 



flow (jLorentzl (jl907lk 



Ladvzhenskaval (|1963f ) ; 



Pozrikidid (| 19921 ) ) . As a reminder, the key 



result of the boundary integral formulation is that a Stokes velocity field u at a point 
x is given by an integral over any closed surface S enclosing x. The surface integral has 
the form 

u(x) = / J(r) ■ [n ■ o-(y)] dS y + f n • K(r) ■ u(y) dS y (3.5) 
is Js 

where y is the point on the surface that you are integrating over and n is an outward- 
pointing surface normal. The tensors J and K arc defined as: 



J(r) 



1/1 rr 



Stt/j \r r 1 



K(r) 



3 rrr 

47r r 5 



r = x - y. 
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Physically, equation (|3.5[) says that the velocity at the interior point x can be expressed 
as a sum over two different kinds of contributions over the the closed surface S. The 
first term on the right hand side of equation (|3.5p corresponds to the component of fluid 
stress normal to the surface S exerted by flows past the enclosing surface. The second 
contribution corresponds to flows into and along the closed surface S (second term on 
the right hand side of equation Q3.5[l ). If the point x lies outside the volume enclosed, 
then the contributions over the closed surface cancel, so that the right hand side of (|3 . 5|) 
sums to 0. A point on the closed surface S is a special case. When the closed surface 
S is continuous and smoothly varying, the velocity at x is simply an average of the 
contribution from S to x if it were in the exterior and the contribution from S to x if 
it were an interior point. As a result, the velocity for a point x on the surface S can be 
written as 

iu(x) = f J(r) • [n • *{y)] dS y + f n • K(v) ■ u(y) dS y . (3.6) 

<j s j s 

These results allow us to derive an integral equation for the time-evolution of the 
interface in our idealized withdrawal problem. To begin, we note that we can enclose 
the entire upper layer as follows: first we define the surface Soo, which has the shape of 
a hemispherical shell with radius Roo- As Roo goes to infinity, the surface 6*00 and the 
liquid interface Si encloses all the liquid in the upper layer. From (|3.5p . we can write the 
disturbance velocity Ui in the upper layer as 

ui(x)= / J (n-o-i) dS y + f n-JT-uidSy. (3.7) 

JSj+S^ JS1+S00 

Since Ui decays rapidly to at infinity, the contribution from the shell Soo approaches 
as Roo is taken to infinity. As a result, the only remaining contributions are due to the 
liquid interface Si, so that we can write the disturbance velocity Ui in the upper layer 
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S R , held at constant p o 

Figure 2. Simplified numerical model of selective withdrawal. The upper and lower liquid layer 
is separated by an interface Si, constrained so that the deflection is nonzero only within a radius 
a. At a finite withdrawal flux, the surface Sr lies entirely within the lower layer. 



MS 

ui(x) = / J (n-o-i) dS y + / n-K-uidSy (3.8) 

JSz JSi 

From (|3.6p . the velocity at a point on the interface is given by 

iu(x) = / J (n-o-i) dS y + / n-Zf-UidSy for x 6 Si. (3.9) 

2 J Si J Si 

Finally, consider a point x on the surface Sr given by z — 0, as sketched in figure O 
At a finite withdrawal flux, the interface is deflected away from the z — plane. As a 
result, points on Sr lie entirely inside the lower layer and therefore are not enclosed by 
the surface Si + Soc , therefore the contribution from the surface integral over Si + ^oo 
vanish at a point x on 5^, or 

= / J ■ (n • <Ti) dS y + / n - K ■ \ii dS y for x in lower layer. (3.10) 

J Si JSi 

Next we consider the volume of liquid in the lower layer enclosed by the closed surface 

comprised of Sr and Si. Again, starting with (|3.6|) . we can write the velocity at a point x 
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on either the liquid interface Si or the surface Sr, henceforth referred to as the "reservoir 
surface", as 

iu(x)=/ J (n-<r 2 ) dS y + f n-K-u 2 dS y + / J • (n • <x 2 ) dS y (3.11) 

2 JSi J Si JSr 

+ n • K ■ u 2 dS y for x e Si or S R 

where n again is an outward pointing surface normal. As a result, for a point x on the 
liquid interface Si , the surface integral over Si in 1|3.11[) has exactly the opposite sign as 
the surface integral over Si in (|3.9[) . Since the velocity is continuous across the interface, 
the two surface integrals involving the K tensor cancel exactly when the two equations, 
(|3.1ip and (|3.9p . for the velocity at a point x on the interface are added together. As a 
result, the velocity on the liquid interface can be re- written as 

u(x) = I J [n-<x]±dS y + I J-(n-cr 2 )dSy (3.12) 

+ / n • K ■ u 2 dS y for x G Si 

Jsr 

where [n • cr]l = n • <j\ — n • <r 2 denotes the jump in the normal stress across Si due to 
surface tension and is equal to 27KI1 where k is the mean surface curvature and n points 
outwards from Si and Sr. Similarly, adding (|3.10[) and (|3.11[) yields an expression for 
the velocity at a point x on the reservoir surface Sr 

iu 2 (x) = / J-[n-cr]ldSy (3.13) 
+ J ■ (n • er 2 ) dS y + / n K ■ u 2 dS y for x e Sr. 

" Sr " Sr 

Equations (|3.12p and (|3.13p together provide an expression for the velocity on the liquid 

interface as a function of the interface shape, and n • <r\s R , the normal stress exerted 

by flow in the lower layer on the reservoir surface Sr. In a full calculation, n • cr\s R 

would need to be obtained separately and then used in (|3 . 1 2[) and (|3 . 1 3(1 to calculate 

the velocity on the liquid interface. Here, we employ a drastic simplification and simply 
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prescribe the normal stress distribution on Sr. Specifically we require that the normal 
stress is a spatially uniform pressure of size po , so that 

n tr 2 |s a =Pon . (3.14) 

This choice of the stress condition (|3.14p is motivated by the observation that, since the 
lower layer is taken to be infinitely deep, the stress distribution far from the interface 
must decay smoothly onto one that corresponds to a stagnant layer. This is simply 
a spatially uniform pressure field. By imposing this distribution at Sr, instead of a 
boundary condition far from the interface, we are essentially approximating the smooth 
decay that would be obtained in a deep lower layer as an abrupt cut-off at Sr. This is 
in the same spirit as approximating the effect of hydrostatic pressure by introducing a 
pinning length a for the interface and turns out be just as effective in capturing the key 
features of the interface evolution. A more realistic simulation which includes a lower 
layer of finite depth, instead of imposing (|3. 14)) at Sr yields essentially the same results 
(section 4.2). 

We are now ready to specify an evolution equation for the liquid interface. The distur- 
bance velocity on the interface is given by 

u(x) = f J ■ n[72«] dS y (3.15) 

J Si 

+ J ■ n p dS y + / n K ■ u 2 dS y for x G Si 

J Sr '-'Sr 

while the velocity on Sr needed to evaluated the third surface integral in (|3.15[) is given 
by (|3.13|) . The total velocity on the interface is a sum of the disturbance velocity and 
the imposed withdrawal flow. Therefore, to update the interface, we use the kinematic 
condition 

dx 

— =u(x)+u cxt forxeSj. (3.16) 
at 

To find a steady-state hump profile for a given withdrawal flux, we start with an initial 
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guess and solve (|3.13[) . (|3.15() . and (|3 . in succession until either a steady-state is 
obtained or the interface deformation grows so large that the interface reaches the sink. 

Finally we nondimcnsionalize the withdrawal problem via the following characteristic 
length-, velocity- and stress-scales 

i = a u=- p=- (3.17) 

(j, a 

The dimcnsionlcss withdrawal flux is then 

Q = ^— (3.18) 

7 

and is equivalent to a capillary number. 

3.2. Numerical Implementation 

We use a C++ code to solve the governing integral equations (|3.13[) . (|3.15p . and (|3.16[) . 
Since the problem is axisymmetric, all the surfaces involved can be represented as effec- 
tively one-dimensional objects by their r and z coordinates in a cylindrical coordinate 
system. The liquid interface is parameterized by mesh points Tj (sj) ,Zj (sj), where Sj 
denotes the arclength along the surface measured from the tip. The interface between 
mesh points is approximated by a cubic spline. This way, the curvature calculated for 
the Laplace pressure term varies continuously between mesh points. The distribution 
and number of points is adapted to the geometric shape of the surface. The algorithm we 
implemented chooses the density of points to be proportional to the mean curvature (in 
most runs As = 0.05/re, where As is the arclength between two adjacent mesh points). 
In addition we require that adjacent mesh points cannot be separated by As larger than 
a maximum value. For most runs, this value was chosen to be of the order of 1/20 in 
dimensionless units. Since the interface at Q = is nearly flat, it is represented by 20 grid 
points. During a typical run, the number of mesh points on the interface is dynamically 

20 



increased to about 60 in order to resolve a dimcnsionless tip curvature of approximately 
75. Doubling the mesh resolutions yields no significant change in the results. 

The constant pressure reservoir surface, Sr, in the setup are represented by a uni- 
form grid of points. The velocity between points are interpolated linearly. The constant 
pressure reservoir surface in the setup is represented by a mesh of 30 uniformly spaced 
points. For numerical reasons, it is simpler to work with a small but non-zero pressure. 
For most runs shown in this paper, po was chosen to be 0.01. As described in section 4.2, 
changing pq results in little relative changes in the results. 

The azimuthal component of the boundary integral equations is in tegrated ana l ytical ly 
and the resulting elliptic integrals are evaluated numerically (see iLee &: Leall (jl982r i). 
The boundary integrals are evaluated using Gaussian quadratures. Singular parts of the 
integral are dealt with separately using finer subdivisions of the integration interval, and 
the resulting algebraic equations for the unknown velocities and stresses are solved by 
LU (lower and upper triangular matrix) decomposition. 

The mesh points on the interface are advanced in time via explicit forward Euler time 
stepping. To ensure that the mesh points remain evenly distributed, only the normal 
component of the velocity is used to update the interface. The length of each timestep 
is chosen to be proportional to the smallest mesh spacing. The validity of this approach 
was tested by doubling the time resolution for a given run and comparing to the original 
results. Generally, results obtained this way differed by less than 1%. We also exper- 
imented with an implicit backward Euler scheme, but the evaluation of the necessary 
Jacobian was prohibitively slow, and the results obtained by the two methods did not 
differ significantly. 

The scheme for finding steady state shapes and critical flow rates follows the experi- 
mental procedure: starting with a static configuration at zero flow rate, Q is increased in 
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small steps and for each step the interface is allowed to relax into a steady state shape. 
Steady state is detected by observing that the modulus of the highest normal velocity on 
the interface, u max n, falls below a certain threshold (10~ 3 in dimcnsionlcss velocity units 
in most cases; lower thresholds were tested without significantly changing the results). 
When it max increases monotonically by a certain factor, we deduce that a steady-state 
hump shape has ceased to exist. In this case, the interval between the last steady state 
flow rate and the current flow rate is divided in half, creating nesting intervals containing 
the critical flow rate Q c . This procedure stops when the size of the interval between Q 
values has decreased below the desired accuracy. To resolve the steady-state evolution 
near Q c , we also choose Q values between those chosen from the bisection and solve for 
the steady-state interface shape at these values. 

Since the successive increases in the imposed withdrawal flux Q becomes very small 
near Q c and the initial velocity on the interface is proportional to the increment in Q, 
the velocity on the interface becomes small near Q c regardless of whether the interface 
shape is stable. We therefore allow the interface to evolve for a longer interval in time 
before checking whether the w max exceeds the threshold value when Q is close to Q c . We 
have also compared our results against runs done with an even longer interval of waiting 
time to ensure that the simulation has converged onto the steady-state solution. 

4. Results 

Using the equations governing the model withdrawal problem described in section 3, we 
calculate numerical solutions for the interface evolution near the transition. Section 4.1 
describes how the curvature k and the hump height h scale with the withdrawal flux Q 
near the transition. This section also addresses the coupling between h and k. Section 4.2 
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shows these results are robust and insensitive to variations in the withdrawal conditions. 
Finally, in section 4.3, the numerical results are compared with those in the experiments. 

4.1. Interface Evolution & Transition 

Time-dependent simulations of equations (|3.13[) . (|3.15[) . and (|3.16[) show that steady-state 
solutions corresponding to hump shapes are found only below a threshold flow rate Q c . 
Above Qc no hump solutions are found. Figure [3] shows a sequence of the calculated 
steady-state shapes. At Q = 0, the static interface is a spherical cap shape determined 
by a balance of surface tension and reservoir pressure po = 0.2. When the imposed 
withdrawal flux Q is small, the interface is weakly perturbed from the spherical cap 
shape. At larger Q, a broad hump develops on the interface. As Q approaches Q c the 
hump height increases slightly but the hump curvature increases dramatically. Above Q c 
the interface develops a finger-like structure which lengthens over time and eventually 
reaches the sink (figure 3 inset). For all Q, the shape of the hump at its tip remains 
smooth, well-fit by a spherical cap with mean curvature k. We identify the disapparance 
of a steady-state hump solution at Q c in the numerical solutions with the hump-to-spout 
transition observed in the experiment. 

In the rest of this subsection, we will analyze one typical set of results, obtained 
with dimcnsionless sink height S = 0.2 and reservoir pressure po set to 0.01. Hump 
solutions obtained for different dimcnsionless sink height S = S/a, po values, and under 
different realizations of two-layer selective withdrawal can be found in section 4.2. We 
can quantitatively characterize the interface evolution by plotting the hump height h 
and the mean curvature k at the tip of the hump as functions of the withdrawal flux 
Q (figure 2]). As Q approaches Q c , the hump height saturates at h c (figure 4a) and the 
curvature saturates at k c . (figure 4b). The insets show that h and n approach their 
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Po=0-2 

Figure 3. Steady-state interface shapes (solid lines) at different Q. For Q — 0, the interface 
assumes a spherical cap shape for which the Laplace pressure balances the reservoir pressure 
Po. As Q increases, the steady-state shape develops a hump towards the sink. The inset shows 
that, above the transition withdrawal flux Q c , the unsteady interface (dashed line) develops a 
finger reaching towards the sink. Here the dimensionless sink height S = S /a = 0.2 and the 
dimensionless reservoir pressure jump po is 0.2. 

saturation values as 

^oc^ ^ Kv ^ 5^%®-. (4.1) 

fl c (eje 

To obtain the best power-law fit, the h c and Q c values were adjusted slightly. The best 
fit is obtained when the values of h c and Q c were increased from the h and Q values 
corresponding to the final steady-state hump shape calculated by 0.03% and 4 • 10 _5 % 
respectively. The same Q c value is used in the k c — k plot and the h c — h plot. The 
k c value is adjusted by less than 0.001% from the final hump curvature value obtained 
numerically. All the adjustments are within the error bars of the simulation. 

The observed square-root dependence of h c — h and k c — k suggests that the transi- 
tion from selective withdrawal to viscous entrainmcnt occurs via a collision at Q c of two 
steady-state hump solutions, one stable and one unstable. In other words, the transition 

from selective withdrawal to viscous withdrawal has the structure of a saddle-node bi- 
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(Q c -Q)/Q c 

Figure 4. Evolution of the calculated steady-state hump shape with withdrawal flux Q. 
(a) Hump height h vs. Q c — Q/Q c . The inset shows (h c — h)/h c approaches saturation as 
V {Qc — Q) /Qc- (b) Mean curvature k at hump tip versus Q. The inset shows (k c — k)/k c 
approaches saturation as \J (Q c — Q)/Q c - The dashed line depicts a square root powerlaw. 



25 



unstable hump 




Figure 5. (a) Stable and unstable hump solutions in steady-state selective withdrawal, (b) 
Saddle-node bifurcation diagram illustrating the evolution of the hump curvature in the selective 
withdrawal regime. 



furcation. Typically, saddle-node bifurcations are associated with solutions that remain 
smooth at the bifurcation point. Therefore the observed square-root scaling suggests that 
the final evolution of the interface as the transition is approached is not organized by an 
approach towards a singular steady-state shape. 

The saddle-node structure of the transition can be rationalized by considering the flow 
in the upper layer converges along the centerline and therefore speeds up away from 
the tip of the hump. Since the hump solution is realized in a time-dependent simulation 
and also observed in the experiment, it is linearly stable. We therefore expect that small 
upwards perturbations of the hump tip simply decay downwards as the hump shape 
relaxes towards the steady-state solution over time. However, if the upwards perturbation 
is very large, so that the perturbed hump tip now lies very close to the sink, then surface 
tension effects will be too weak to pull the perturbed interface downwards. Instead, the 
hump tip will be drawn into the sink. This suggests there exists an intermediate critical 
shape perturbation that neither decays nor grows upwards but simply remains steady 
over time, as illustrated in figure [5^,. This critical shape perturbation would correspond 
to an unstable hump solution. 
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Figure 5b illustrates the evolution of the hump curvature as a function of the flow 
rate for stable and unstable hump solutions. At low withdrawal flux values, only a very 
large shape perturbation can cause the interface to be drawn into the sink, we therefore 
expect the stable and the unstable hump shapes to be widely separated. This means 
the unstable hump must lie close to the sink and, as a result, has a very curved tip. 
At moderate withdrawal flux values, the perturbation size required to destabilize the 
hump solution becomes smaller, implying that the unstable and stable solutions now lie 
closer to each other. Concurrently, the curvature of the stable solution increases, and the 
curvature of the unstable solution decreases. Near transition, even a small perturbation 
causes the interface to become unstable and grow towards the sink. This suggests that 
the two solutions lie very close each other and are nearly identical. Therefore it is very 
reasonable to expect the two solutions to become identical, or coincide, at Q c , thereby 
bringing about a saddle-node bifurcation of the hump solution. The square-root scaling 
corresponds to a smooth merging of the two solutions, so that the k(Q) curve at Q c has 
a shape of a parabolic lying on its side. 

The scaling dynamics associated with how the hump solution saturates as Q approaches 

Q c is consistent with a saddle-node bifurcation at Q c . This is a generic mechanism for 

transition for one type of solution from another in a nonlinear problem. However, the 

quantitative analysis in figure [4] shows one surprising and atypical feature: the hump 

height and curvature saturate towards the final scaling behavior at very different 5q 

values. Specifically, the square-root scaling in h c — h is evident throughout its evolution, 

even at large Sq. In contrast, the square- root scaling for k c — n becomes evident only 

when Sq has decreased below 10 -3 . This behavior suggests that the hump shape does not 

approach the transition uniformly. To test this idea we plot the hump radius at z = h/2 

as a function of Q in figure [6] We find that the radius at the half-height saturates to the 
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Figure 6. Evolution of the hump radius at z = h/2 as a function of the withdrawal flux Q. 
The inset shows how the radius at half height reaches its saturation value. 

final scaling behavior later than h, but earlier than n as Sq approaches 0. This behavior 
indicates that, as Q approaches Q c , the overall shape of the hump, e.g. the hump height 
or its lateral extent, saturates first, followed by features on smaller lengthscales. The 
shape of the hump at its tip, which corresponds to a feature on the smallest lengthscale, 
saturates last. This cascade of events is more typical of an approach towards a singular 
shape, in which the evolution of features on different lengthscales are nearly decoupled, 
so that features on smaller lengthscales saturate later than features on large lengthscales. 

To get some insight into this unusal hybrid character of the interface evolution near 

transition, we plot n as a function of h in figure [7] At Q = 0, the interface shape is 

given by a balance of Laplace pressure and reservoir pressure pq. For po — 0.01, the 

shape is a nearly flat spherical cap. When the interface is only weakly perturbed from 
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the Q = shape, both k and h increase linearly with Q, and therefore k oc h. This linear 
regime persists until h at 0.02. Beyond this point, k increases much more rapidly than h. 
Remarkably, this large increase of k relative to the increase in h is well-approximated by 
a simple mathematical expression, as evident from the inset for figure [7] In the nonlinear 
deflection regime, 

1h(k) cx h . (4.2) 
This logarithmic coupling indicates that in the nonlinear regime, the hump solution 
evolves so that small increases in hump height corresponds to large increases in the hump 
curvature. This does not corresp o nd to a nearly continuous transition, as suggested by 



Cohen & Nagel (jCohen fc Nagell (|2002h ). since the hump curvature k never decouples 
from the hump height. Instead, our results show that k remains weakly coupled to h and 
does not diverge as the tranisiton is approached. To gain some insights into the origin 
of this unusual coupling, we next investigate how changes in the withdrawal conditions, 
in particular the sink height S and the reservoir pressure po, affect this relation between 
the hump curvature and the hump height. 

4.2. Interface Evolution under Different Withdrawal Conditions 

Naively one may expect that a logarithmic coupling between the hump height and the 
hump curvature would be easily perturbed by changes in the boundary conditions, since 
logarithmic coupling is often observed at the transition between different kinds of power- 
law scaling as some system parameter is varied. In addition the withdrawal process 
analyzed numerically is highly idealized and may therefore give results which are not 
generic. To address these concerns, we analyze numerical results for the steady-state 
interface shapes under a variety of withdrawal conditions. Surprisingly, we find the log- 
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Figure 7. Mean curvature at the hump tip versus the hump height. Inset shows the same data 
on a semi-log plot. When the deflection is sumcently large, the hump height h has a linear 
dependence on ln(ft). 

arithmic coupling (|4.2[) to be a robust feature of equal-viscosity selective withdrawal in 
all the numerical solutions. 

To determine the effect of the reservoir pressure po on the transition structure we 
analyzed the numerical solutions for values of po ranging from —0.3 to 0.3. Larger po 
values were not investigated because they result in static shapes whose height is greater 
than the sink height even at Q = 0. In figure [S] we plot k c and h c , the hump curvature 
and hump height at transition versus pq. We find that the precise value of po makes 
little difference to the values of re c and h c . Also, all the runs show the same logarithmic 
dependence of k on h near the transition. 

In contrast, figure [9] shows that varying the sink height S strongly affects k c and h c . 

As the dimensionlcss sink height S becomes much smaller than 1, or S <C a, the hump 

curvature at transition, n c , increases as 1/S while h c decreases as S. In the opposite limit, 

as S becomes much larger than 1, or S > o, both n c and h c approach constant values. 
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Figure 8. Calculated hump curvature k c and hump height h c at transition as a function of 
the reservoir pressure po. The sink height 5=1. 

These trends can be understood by analyzing the limiting cases where the sink is either 

very far or very near to the interface. When the sink is placed very close to the interface, 

the effect of the pinning at a becomes irrelevant since the interface is already flat at 

distances O(S) where the flows become very weak. Consequently, h c and n only depend 

on S when S <C 1. In the opposite limit, where S ^> a, the flow at the interface becomes 

almost uniform because the lengthscale over which the imposed withdrawal flow varies is 

much larger than a. Under these conditions the interface shape is primarily determined 

by the pinning condition at a and is insensitive to changes in the sink height. 

Figure [9] also shows that the variation in k c tracks the variation in h c almost perfectly 

as S is varied. Consequently, even though driving the withdrawal with a sink closer to the 

interface produces a larger hump curvature at transition, it does not prduce a sharper 

hump tip. This is because the hump height, and similarly the overall deformation at 

31 




1000 




Sink height S 

Figure 9. Calculated hump curvature k c and hump height h c at transition different sink heights 
S. The reservoir pressure po = 0.01. The dashed line in the k c vs. S plot shows a 1/S dependence. 
The dashed line in the h c vs. S plot shows a S dependence. 

transition, is proportionally smaller as well, so that the product n c ■ h c characterizing the 
relative separation of lengthscale between the hump tip and the hump height remains 
roughly constant with 5*. Changing the sink height primarily changes the absolute length- 
scale of the steady-state deformation at the onset of entrainment. It does not change the 
shape of the deformation significantly. 

This observation suggests that dividing all lengths by h c should scale out most of the 
variation observed at different sink heights. In figure [TTH wc plot h c -n versus h/h c obtained 
from numerical solutions where S = 0.05, 0.2, 1, 10 and 50 representing a 10 3 variation 
in the sink height. We find that the curves are indeed brought close together by this 
rescaling. For example, while the actual value of k c changes by a factor of 100 between 

S = 0.05 and S = 50, the rescaled value h c ■ k c changes by less than a factor of 4. We 
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Figure 10. Calculated hump curvature k versus hump height h at S = 0.05, 0.2, 1, 10 and 50. 

do however, observe small changes in the slope and intercept of the logarithmic curves 

even when k and h are rescaled. Nevertheless, from these results, we conclude that the 

logarithmic coupling is a robust feature of our idealized model of two-layer withdrawal. 

Changing either the boundary condition parameter p , or the forcing parameter S within 

the simple model of withdrawl produces no qualitative change in the evolution of the 

steady-state hump shape as a function of the withdrawal flux. 

Given this robustness with respect to variation in system parameters in the simple 

model, we next analyze the steady-state interface obtained under different realizations 

of selective withdrawal. We find that even these different realizations, which correspond 

to different boundary conditions and/or imposed withdrawal flows, fail to produce a 

qualitative change in how the steady-state hump shape evolves. Figure [TT] shows the 

interface shapes obtained in two different sets of calculations, each with a lower layer 

of finite depth, instead of the infinite depth assumed in the model calculation. In the 

first case, the lower layer is contained in a cylindrical cell of radius a\ and depth a\ 
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Figure 11. (a) Steady-state interface shapes obtained in selective withdrawal from a lower 
layer of finite depth. The cell depth equals the cell radius a 1 and the sink height S/a 1 = 0.5. (b) 
Steady-state interface shapes obtained in selective withdrawal from a thin lower-layer overlying 
a solid plate of radius 0,2 . The sink height S/a2 = 0.2. 

(figure II lb .) . This corresponds to a situation where the lower-layer is confined equally 
in both the radial and vertical direction. The toroidal recirculation established in the 
lower-layer is therefore much smaller in extent, thereby resulting in an O(l) change in 
the viscous stresses exerted on the interface by the flow in the lower layer. 

For this realization, the boundary integral formulation uses a closed surface comprised 
of the liquid interface 57, the side walls of the container S S id c and the bottom wall of the 
containter Sb- The disturbance velocity u on the interface is given by 

iu(x) = f J n [t2k] dS y + f J n er sido dS y (4.3) 

+ I J ■ n er b dS y for x e Si 

Js b 

where y is the point on the surface that you are integrating over, <r s idc and <Xb correspond 
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to the normal stress exerted by the flow on the side wall and bottom wall of the container. 
The two stresses are found by solving the equation 



This corresponds to enforcing no-slip boundary conditions on the side walls and the 
bottom walls of the container. 

The interface shapes shown are the calculated results and the cell radius ai is used as 
a characteristic lcngthscale in nondimcnsionalizing various quantities. As the withdrawal 
flux Q increases and the hump height increases, liquid is drawn from the lower layer into 
the hump. Since volume in the lower layer is conserved, gathering extra liquid into the 
hump forces a shallow dip in the interface shape to develope some distance away from the 
ccntcrline. Thus, the hump profile no longer flattens monotonically with r, in contrast 
to results obtained assuming an infinitely deep lower layer. However, these qualitative 
changes to the overall shape do not significantly change how the interface evolves with 
Q. The shape transition still corresponds to a saddle- node bifurcation (figure [T2]) . A 
comparison of the rescaled h c ■ k versus h/h c curve obtained for a cell of finite depth 
and the analogous curve obtained assuming an infinitely deep layer in figure [T3] show the 
same trend. 

Figure [Tib shows what happens when we change the condition in the lower layer more 
drastically. In this set of calculation, a thin layer of the lower layer liquid overlies a solid 
plate. The plate radius ai is used as a characteristic lcngthscale in the nondimensional- 
ization. The liquid is assumed to wet the solid substrate so that the plate is covered by 
liquid from the lower-layer at all Q. The formation of a hump causes most of the liquid 
in the layer to drain into the hump, leaving only a very thin layer covering the rest of 




(4.4) 
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Figure 12. Scaling of the calculated hump height h near the transition flow rate for two- layer 
withdrawal. Set (a) corresponds to withdrawal from a cylindrical container. Set (b) corresponds 
to withdrawal from a thin lower layer which is resting on a solid plate. 



the solid plate. In the selective withdrawal regime, the withdrawal flow requires liquid 
in the lower layer to recirculate, moving along the interface towards the hump and then 
along the lower solid surface away from the hump. When the lower layer is thin, the vis- 
cous resistance associated with the steady-state recirculation is amplified dramatically. 
In this realization, therefore, the stress-balance on the steady-state interface is strongly 
perturbed from that obtained for withdrawal with an infinitely deep lower layer. 

For two-layer withdrawal with a thin layer of the lower liquid, the closed surface used 
in the boundary integral formulation consists of the liquid interface Si and the bottom 
wall of the containtcr Sb- The equation for the velocity on the interface has the same 
form as (|4.3j) except that the terms associated with the side wall surface SWde are absent. 



30 



The normal stress on the solid surface Ob satisfies equation (|4.4[) . without the sidcwall 
contribution. 

Comparing against results from the simple model and results for withdrawal from a 
finite container, we can see that reducing the thickness of the lower layer has an effect 
on the steady-state shape of the interface obtained. The hump shape is significantly 
more conical, and the flattening of the interface at large r results primarily from volume 
conservation, rather than pinning at 02 or the decay of the imposed withdrawal flow 
away from the sink. However, as evident from figure [12] and figure 1131 the qualitative 
features are unchanged. The transition still corresponds to a saddle-node bifurcation. 
The evolution of the hump height h retains a logarithmic coupling to the hump curvature 
k. The main difference between withdrawal from an infinitely deep layer and withdrawal 
from a very thin layer overlying a solid plate is a shift in the slope of the rescaled k versus 
h curve. 

We have also conducted a series of calculations which include the effect of hydrostatic 
pressure explicitly. This requires that we change the form of the normal stress jump 
[n • a]_ across the liquid interface Si in the J-integral of equations (|3.13[) and (|3 . 1 5[) to 

[ 7 2k - Apg • y]n (4.5) 

where Ap is the density difference between the two layers, g the gravitational acceleration 
and y the location of a point on the liquid interface. Choosing the capillary lcngthscale 
£ 7 to be less than the pinning radius a allows the interface to flatten out at large radial 
distances due to stratification, as occurs in the experiment, instead of due to the decay 
of the withdrawal flow or due to surface tenson effect. Results for £ 7 /a = 0.1 and 0.2 
show that introducing stratification explicitly results in slight qualitative changes in the 
hump shape but the same trend in the rescaled k versus h curve ( figure fT3|) . Finally, to 
assess the influence of details of the geometay of the withdrawal flow on the logarithmic 
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Figure 13. Curvature versus hump height under different withdrawal conditions. Open circles 
correspond to withdrawal with a lower layer of finite depth contained in a cylindrical cell. The 
cell depth equals the cell radius eti and the sink height S/a\ = 0.5. The crosses correspond to 
withdrawal of a thin layer of liquid overlying a solid plate of radius and sink height S/a,2 = 0.2. 
The open triangles correspond to withdrawal from a stratified layer. The capillary lengthscale 
l-y/a = 0.1 and the rescaled sink height is 0.2. The solid diamonds correspond to withdrawal 
driven by a vortex ring of radius a u = a, dimensionless strength u and height S/a = 0.2. The 
K-h curve obtained for S = 0.2 and po = 0.01 (solid line) is also shown for comparison. 



coupling (|4.2[) . wc changed the imposed withdrawal flow from a sink flow (|3.2[) , to that 

associated with a vortex ring of size a u = a and strength loq. The axisymmetric velocity 

field is given in terms of a stream function <f> 

1 50 1 90 

u r = -— u z = --— (4.6) 
r or r oz 

and the stream function has the form 

Iff f 27T cos 9 AO 



<t>(r,z) = -7- rr'u;(r',z'Wdr' / ==__=_ (4.7) 

47T J J J Q y/(z - z ') 2 + r 2 + r' 2 - 2rr' cos 9 

where the vorticity distribution u>(r' , z') has the form LOoS(a^ — r)8(S — z') for a vortex 
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ring. This is a better approximation of the withdrawal flow imposed in the experiment, 
both because there is no net volume withdrawn and because the finite size of the vortex 
ring mimics the finite diameter of the tube. The rescaled k versus h curve obtained from 
numerical solutions for withdrawal driven by a vortex ring again has the same logarithmic 
coupling between the hump height and the hump curvature near the transition. 

In all the different realizations of selective withdrawal, the evolution of the steady-state 
interface approaches the form (|4.2|) for large deformations. As a result, h c ■ k c , the relative 
separation of lcngthscales characterizing the hump shape at transition, changes little 
even when the forcing and/or the boundary conditions are changed drastically. In other 
words, regardless of what type of flow we use to drive the transition, the degree of viscous 
resistance from flow within the lower layer, or whether the interface flattens out on large 
lengthscale due to stratification or surface tension effects, the hump shape at transition 
never becomes significantly sharper than that first calculated in the simple model. We 
conclude from this that the interface shape at Q c produced in viscous withdrawal with 
two liquid layers of equal viscosity cannot be brought closer to a singular shape via 
changes in the boundary conditions. This contrasts sharply with theoretical results on 
viscous entrainment when the entrained liquid is far less viscous than the exterior, which 
show that interface shap e at Q c can be made singular via changes in the boundary 
conditions (jZhanj (|2004f )). 

Finally, the robustness of the logarithmic coupling even under drastic changes in the 
withdrawal conditions suggests that (|4.2p in fact corresponds to a generic feature of selec- 
tive withdrawal from two layers of equal viscosity. This leads us to re-examine experimen- 
tal measurements of k and h, previously interpretted in terms of a continous evolution 
towards a change in interface topology. In the next subsection we compare experimental 

results against numerical results obtained using the simplest model of withdrawal and 
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show that the logarithmic coupling (|4.2[) is also a good description of the steady-state 
interface evolution in the experiments. 

4.3. Comparison with Experiment 
Here we compare three key r esults from the numeri cs against the measurement previously 



obtained by Cohen & Nagel (jCohen fc Nagell (|2002f )). First, we analyze the h(Q) obtained 



in the experiment and show that the measurements arc consistent with the interpretation 



that h c — h scales with \/Q c — Q near the transition. Measurements of the hump curvature 
also agree well with the calculated values and begin to saturate as Q c is approached. 
Finally we compare the measured n versus h curves against the calculation and show 
they also agree. 

To make the comparison, we chose measurements from five different experiments, span- 
ning the full range of tube heights used. Figure [TJ] shows how the measured hump height 
saturates as Q approaches Q c . As was done with the numerical results, in generating 
figure [14] we allowed ourselves to vary h c and Q c values within the experimental error 
bars, which are about 5%, in order to generate the best power-law fits for the measure- 
ments. For four sets of the data, the saturation behavior is completely consistent with a 
square-root scaling with 5q. The set with the largest tube height (S p = 0.830 cm) shows 
a slight difference in the scaling behavior. We do not fully understand the origin of this 
discrepancy. We speculate that withdrawal experiments at larger tube heights require 
running the pump at higher flow rates. This may create more noise and/or signficant in- 
ertial effects, hence resulting in a discrepancy with experiments performed at lower tube 
heights. Overall the agreement shows that the hump height in the experiment experiences 
a saddle-node bifurcation at Q c . 

Next we compare measurements of the hump curvature against calculated values. Since 
the hump curvature satures at a much smaller value of 5q, below the dynamic range 
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Figure 14. Rescaled hump heights (h c — h)/h c versus rescaled withdrawal flux (Q c — Q)/Q c 
for measurements obtained with 5 different tube heights. The tube heights are S p = 0.260 cm, 
0.381 cm, 0.518 cm, 0.712 cm and 0.830 cm. The capillary lengthscale in the experiment is 0.3 
cm and the layer viscosity ratio (lower/upper) is 0.86. The calculation results are for S/a = 0.2 
and reservoir pressure p = 0.01. 



of the experiment, we cannot compare the saturation dynamics, which is independent 
of the details of the numerical model, directly. Instead we compare the ii(dq) curves 
obtained experimentally against the results obtained using our minimal numerical model, 
described in section 3.1, with S = S /a — 0.2. This is done in figure [T5l The Q c values 
which produced the best power-law fit for the experimental data in figure [H] are used in 
figure IT5l (More infor mation on the differenc e between our analyses and those performed 



in the original papers (jCohen fc Nage 



( 20021 )) can be found in the appendix) As was done 



for numerical results obtained with different sink heights, we account for the change in the 
absolute size of the hump at different tube heights by rescaling k by h c , the hump height 



41 




10 



J3 



10" 



10" 



10" 




simulation 
c^S=0.255cm 
^S=0.381cm 
<^S=0.5 175cm 
^S=0.7 124cm 
>~< S=0.8295cm 



10" 



10 10"' 

(Q-Q)/Q, 



10" 



10 



Figure 15. Tip curvature versus rescaled flow rate (Q c — Q) /Q c for the numerical model and 
experimental data. The tip curvature has been non-dimensionalized by h c . Range of the numer- 
ical results have been truncated to display the comparison clearly. 



at transition. This causes the different curves associated with different tube heights to 
collapse onto roughly a single curve. Significantly the collapsed curve shows some evidence 
of saturation as 5q approaches 0. In order to display the comparison clearly, the range of 
value for both k ■ h c and (Q c — Q)/Q c have been truncated. The calculated curve goes 
through the experimental values and shows exactly the same trend. Note that our choice 
of S = 0.2 for the numerical results is primarily a matter of simplicity, since that is the 
set of results discussed in section 4.1. From section 4.2 we have seen that k c ■ h c varies 
only weakly with S, so we could have used any 0(1) value for the dimcnsionaless sink 
height S and gotten good agreement in this rescaled curvature plot. We conclude from 
the good agreement that the hump curvature also saturates in the experiment. 

Finally we plot the rescaled measured hump curvature as a function of the rescaled 
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Figure 16. Comparison between measurements at five different pipette heights and calculation. 
Note the range of k ■ h c has again been truncated to display the overlap with measurements 
more clearly. 



hump height h/h c (Figure [T6|) . The measurements at the largest tube height (S p = 0.830 
cm) are displaced from the four other sets. Also, the measured curves for the three, larger 
tube heights show a slight but systematic upturn at large h/h Cl when compared against 
the calculated results. The two data sets for the lower tube heights appear to follow a 
logarithmic relation. While it is unclear what causes these discrepancies, overall there is 
good agreement between the calculated curve and the measured curves. The agreement 
shows that the simple model used here captures the main features in the evolution of the 
steady-state interface shape. 

Previously, the large increase in the hump curvature near transition was interpret- 
ted in terms of the steady-st ate interface app r oaching a stea d y-stat e singularity as the 



withdrawal flux is increased (jCohen fc Nage 



(120021 ): 



Cohenl (J200J)). The good agree- 
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merit between the numerical solutions and the experimental measurements show that 
this large increase in hump curvature can also be understood to result from a weak cou- 
pling between the overall hump shape, characterized by the hump height, and the shape 
of the hump tip, characterized by the hump curvature. The simplest way to account for 
all the results, numerical and experimental, obtained for equal- viscosity withdrawal is to 
assume that no singular hump solution exists for equal- viscosity withdrawal and that the 
logarithmic coupling (|4.2p is a generic feature of selective withdrawal. The absence of a 
singular solution for equal-viscosity withdrawal would explain why even drastic changes 
in the withdrawal conditions failed to introduce qualitative changes in the k versus h 
curves calculated, and why the product k c ■ h c characterizing the sharpness of the hump 
profile always remains below 20. The saddle node structure of the transition and the 
logarithmic coupling between h and k also account in a natural way for the two features 
which appear puzzling in the context of a nearly continuous shape transition: the appar- 
ent transition cut-off in the hump evolution and the large changes in the hump curvature 
relative to changes in hump height. 

While our results simplify the view of the steady-state interface evolution significantly, 
they do not offer a full explanation. We do not know why the logarithmic coupling should 
arise in selective withdrawal in the first place. Section 5 offers a brief, and speculative, 
discussion of its possible origin, focusing in particular on a comparison between previously 
proposed singular solutions and the logarithmic relation found in our numerical solution. 

5. Discussion 

A natural way to reconcile the logarithmic evolution of the hump height relative to 
hump curvature (|4.2[) obtained here with the previously proposed idea that the sharpen- 
ing of the steady-state hump tip is associated with an approach towards a steady-state 
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singularity is to suppose that the interface evolution is indeed organized by the existence 
of a steady-state singularity, as suggested by previous studies. However, the singular 
shape is not nearby in the sense that it is attained at Q just above Q C7 where a transi- 
tion from selective withdrawal to viscous entrainment takes place. Instead, the singular 
solution corresponds to an isolated singularity, existing at a value Q* much larger than 
Q c . Correspondingly, the height of the singular hump is much larger than h c , the 
hump height at transition. In other words, the steady-state interface evolves towards a 
singular shape at Q» as the withdrawal flux is increased. However, this evolution is cut- 
off at a withdrawal flux Q c far below by a saddle-node bifurcation. While unusual, a 
steady-state cusp solution correspon ding to an isolated singularity h as been found for 2D 



selective withdrawal of inviscid fluid (jTuck fc Vanden-Broeckl (| 1984T ) ) so this is a possible 
scenario for the steady-state evolution. 

Under these assumptions, the power-law scaling 



K 

c 



(5.1) 



- h / 

where C is a characteristic curvature scale, indeed assumes an approximate logarithmic 
form. Since h c -C h*, the ratio h/h* is small over the entire range. Taylor series expansion 
of the natural log of (|5.ip yields 



/ K \ 



l*(^J=-/31n(l--]=/?(- + -(- 



1 / h 



(5.2) 



2 V /ijjt 

When h c /h* is sufficiently small, the higher order terms are negligible, and we obtain a 
logarithmic relation between n and h 



h* 
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kr 



(5.3) 



According to Cohen & Nag el, the singular shape i s characterized by a scaling expo- 



nent P between 0.7 5 and .82 ( Cohen fc Nagell (|2002f 0. In the long- wavelength analysis 



of spout formation (jZhand (|2004f )). Zhang ^gund that the steady-state spout shape ap- 



proaches a singular, conical shape. This would correspond to /3 = 1. More recently, exper- 
imental and scaling analysis by Courrech du Pont and Egger s show that the cusp observed 



in a vi scous drainage experiment approaches a conical shape (jCourrech du Pont fc Eggers 



(|2006l )). with (3=1. Given these values, consistency requires that the coefficient [3h c /h* 
must be much smaller than 1, since (|5.3[) is derived under the assumption that h c /h t <C 1. 
We can check this against the numerical results by fitting the calculated k versus h curves 



with the form 



In 



= b 



(5.4) 



where k is the extrapolated intercept at h = and the slope b = (3h c /h* in (|5.3p . Fits 
to the calculated curves for S = 0.05, 0.2 and 50 yield, respectively, b values 5.5, 5.9, and 
1.9, all above 1. Consequently, h c /h* is O(l) which invalidates the Taylor expansion in 
equation ()5.2() . This clearly rules out the possibility that the logarithmic relation reflects 
the existence of an isolated singular hump solution for equal-viscosity withdrawal. 

A number of future studies would complement this one. It would be useful to explicitly 
calculate the unstable hump sol ution, g i ven it s relation with the stability of the hump so- 



lution s. In light of recent works (jChaiebl ((2004); 



Zhang||200J); 



Courrech du Pont fc Eggers 



(|2006h ) suggesting that the steady-state interface indeed evolves towards a singular shape 
when the lower layer is much less viscous than the upper layer, it would also be worth- 
while to analyze the selective withdrawal regime when the two liquid layers have unqual 
viscosities. We expect that the final transition from selective withdrawal to viscous en- 
trainment occurs via a saddle-node bifurcation even when the liquid layers have unequal 
viscosities. However, the limiting situation where the lower layer has a vanishingly small 
viscosity relative to the upper layer corresponds to a qualitative change in the entrain- 
ment dynamics since only the flow in the upper layer is significant. This then suggests a 

qualitative change in the steady-state shape evolution as well. 
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6. Conclusion 

In conclusion, this paper presents a numerical study of a simple example of flow-driven 
topological transition: the transition from selective withdrawal to viscous cntrainmcnt. 
To uncover the fundamental nature of the topology change, we focus on the selective 
withdrawal regime, and analyze how the steady-state shape of an interface between two 
immiscible and viscous liquid layers changes as the imposed flow strength increases. We 
analyze a model problem in which the interface deformation is solely controlled by vis- 
cous stresses and surface tension. Its results show that the hump solution corresponding 
to the selective withdrawal regime experiences a saddle-node bifurcation at a threshold 
flow rate. Above the threshold, no hump solution exists. As observed in the experiments, 
the hump curvature increases dramatically near transition while the hump height in- 
creases only weakly, but the increase corresponds to a logarithmic coupling between the 
hump height and the hump curvature. This suggests that the evolution of the hump tip 
remains weakly coupled with the evolution of the overall hump shape as the transition is 
approached. Moreover, we find that this logarithmic coupling is robust. Numerical solu- 
tions for different realizations of viscous withdrawal, as well as measurements, all show 
a logarithmic coupling between hump height and hump curvature. 
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Appendix A. 

As was seen in previous sections, compared to the hump height, the tip curvature 
shows evidence of saturation towards a final value only when Q is very close to Q c . 
As a consequence, when the entrainment threshold is approached, the tip curvature can 
appear to diverge while the hump height reaches a saturation value. In particular, when 
k is plotted against (Q c — Q)/Q instead of {Q c — Q)/Q c , it is very difficult to distinguish 
between a continued power law and a turn over into saturation with the available range of 
experimental data. This is because even a slight shift in the value of Q c can straighten out 
the curves. The apparent power law in the {Q c — Q)/Q range between 1 and 10 -2 seems to 
be an artifact of the fact that, when Q is small, the denominator of (Q c — Q) /Q is close to 
Q c and almost constant, thus the entire expression scales roughly like 1/Q. Furthermore, 
for small flow rates Q the system response to the outside flow can be assumed to be 
almost linear, with a linear dependence of k on Q. This would then asymptotically give 
rise to a powerlaw with an exponent of —1. 
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